

In [1]:

    
#import functions
%pylab inline

# from MyUnits import *
from MyFunctions import *
from qutip import *

# from MyQubit import *
# import mpld3
import multiprocessing as mp
import itertools
import datetime









    



Populating the interactive namespace from numpy and matplotlib



In [2]:

    
%%javascript
IPython.load_extensions('usability/codefolding/main');
IPython.load_extensions('toggle_all_line_number.js');



In [3]:

    
import scipy.constants as sc



In [100]:

    
import time
import datetime
from scipy.optimize import fsolve

Standard Hamiltonian

$$ \lambda \left({{b}^\dagger} {\sigma_-} + {b} {\sigma_+} - \frac{g^{2} {{b}^\dagger} {\sigma_-}}{2 \Delta^{2}} - \frac{g^{2} {b} {\sigma_+}}{2 \Delta^{2}} - \frac{g^{2} \left({{c}^\dagger}\right)^{2}}{\Delta^{2}} {\sigma_-} {b} - \frac{g^{2} \left({c}\right)^{2}}{\Delta^{2}} {\sigma_+} {{b}^\dagger}\right) + w_{mr} {{b}^\dagger} {b} + {{c}^\dagger} {c} \left(\omega_{c} + \frac{g^{2} {\sigma_z}}{\Delta} - \frac{g^{2} \lambda}{\Delta^{2}} {{b}^\dagger} {\sigma_-} - \frac{g^{2} \lambda}{\Delta^{2}} {b} {\sigma_+}\right) + {\sigma_z} \left(\frac{w_{q}}{2} + \frac{g^{2}}{2 \Delta} + \frac{g \lambda}{\Delta} {{b}^\dagger} {c} + \frac{g \lambda}{\Delta} {{c}^\dagger} {b}\right) $$



In [147]:

    
# Time Hamiltonian
def Ht(t, args):
    #
    # evaluate the hamiltonian at time t. 
    #
    
    H0 = args['H0']
    c = args['c']
    cDag = args['cDag']
    A = args['A']
    
    w  = args['w']
    
    
    s = args['s']
    sx = args['sx']
    
    

    return H0 + A * (c + cDag + 2 * s * sx)*cos(w*t) #(a * exp(1j*w*t) + aDag * exp(-1j*w*t))

def func(a,g,L,w_c,w_nr,w_q):
    s,alfa = a
    return (s+ g/(w_c+w_q),alfa+ L/(w_nr+w_q))



In [148]:

    
# Calc Spectrum
def calc_spectrum_6(N,M,P, w_c,w_nr, w_q,L,g,A,w=0, **kwargs):
    
    # dispersive Qubit CPW NR 
    Delta = w_q - w_c
    delta = w_q - w_nr
    # qubit operators
    
    sm = tensor(create(2),qeye(M),qeye(P))
    sz = tensor(sigmaz(),qeye(M),qeye(P))
    sx = tensor(sigmax(),qeye(M),qeye(P))
    nq = sm.dag() * sm
    xq = sm + sm.dag()
    I = tensor(qeye(2), qeye(M),qeye(P))
    
    
    # mechanical resonator operators
    
    b = tensor(qeye(2),destroy(M),qeye(P))
    n_b = b.dag() * b
    x_b = b.dag() + b
    p_b = b - b.dag()
    
    
    # CPW operators
    
    c = tensor(qeye(2),qeye(M),destroy(P))
    n_c = c.dag() * c
    x_c = c.dag() + c
    p_c = c - c.dag()
    
    # Identity
    
    I = tensor(qeye(2),qeye(M),qeye(P))
    
    
#      Hamiltonian
    if 'Fred' in kwargs:
#         s = -0.01999 + 0.00393001*w_q - 0.000689326*w_q**2 + 0.0000851242*w_q**3 - 4.92054e-6*w_q**4
#         alfa = -0.000285256 + 0.0000784123 * w_q - 0.0000178779 * w_q **2 + 2.60578e-6* w_q**3 -1.65641e-7*w_q**4

        s,alfa = fsolve(func,(0,0),args=(g,L,w_c,w_nr,w_q))
            
        
        w_qt = w_q * exp(-2* s**2 - 2 * alfa**2)
        
        
        H1 = w_c * c.dag()*c + w_nr * b.dag()*b 
        
        H2 = w_qt/2 * sz *(1 - 4 * s**2 * c.dag() * c)* (1 - 4 * alfa **2 * b.dag() * b )

        H3 = (w_qt/(w_c + w_qt))* 2 * g * ( c.dag()*sm + c*sm.dag()) 
        + (w_qt/(w_nr + w_qt))* 2 * L * ( b.dag()*sm + b*sm.dag()) 
    
    
    
    
    
    
  
    
    
    
    
#     Colapse Operators
    
    c_op_list = []
    
    kappa_n = 0.0005 # cavity
    
    gamma_rel = 0.0001 # qubit
    gamma_dep = 0.002 # qubit
    
    Gamma_m = 0.01 # MR
    
    Ta = 60e-3 #k
    Tb = 60e-3 #k
    Tq = 30e-3 #K
    
    n_th_a = 1/(exp(sc.h*w_q*1e9/(sc.k*Ta)-1))
    
    n_th_q = 1/(exp(sc.h*w_q*1e9/(sc.k*Tq)-1))
    
    if Tb == 0:
         n_th_b = 0
    else:
        
        n_th_b = 1/(exp(sc.h*w_nr*1e9/(sc.k*Tb)-1))
    
    # cavity
    c_op_list = []

    rate = kappa_n * (1 + n_th_a)
    if rate > 0.0:
        c_op_list.append(sqrt(rate) * (c+s*sx))

    rate = kappa_n * n_th_a
    if rate > 0.0:
        c_op_list.append(sqrt(rate) * (c.dag()+s*sx))

    rate = gamma_rel * (1 + n_th_q)
    if rate > 0.0:
        c_op_list.append(sqrt(rate) * (1/2 * (sx 
                                              + exp(-2*s**2 - 2 * alfa**2)*(
                        (1 - 4 * s**2 * c.dag() * c)
                        * (1 - 4 * alfa **2 * b.dag() * b )
                        *(sm-sm.dag())
                        -2*s*(c.dag()-c)*sz - 2*alfa*(b.dag()-b)*sz) 
                                             )))

    rate = gamma_rel * (n_th_q)
    if rate > 0.0:
        c_op_list.append(sqrt(rate) * (1/2 * (sx 
                                              + exp(-2*s**2 - 2 * alfa**2)*(-
                        (1 - 4 * s**2 * c.dag() * c)
                        * (1 - 4 * alfa **2 * b.dag() * b )
                        *(sm-sm.dag())
                        +2*s*(c.dag()-c)*sz + 2*alfa*(b.dag()-b)*sz) 
                                             )))

    rate = gamma_dep / 2 * (1 + n_th_q)
    if rate > 0.0:
        c_op_list.append(sqrt(rate) * (sz*exp(-2*s**2-2*alfa**2)
                                       * (1 - 4 * s**2 * c.dag() * c)
                                       * (1 - 4 * alfa **2 * b.dag() * b )))
    
    rate = Gamma_m * (1 + n_th_b)
    if rate > 0.0:
        c_op_list.append(sqrt(rate) * (b+alfa*sx))

    rate = Gamma_m * n_th_b
    if rate > 0.0:
        c_op_list.append(sqrt(rate) *( b.dag()+alfa*sx))    
        
    
# Solution Type    
#     if 'dispersive' in kwargs:
               
#         H0 = H1 + H2a + H3 + H4 #+ H5
#         rho = steadystate(H0,c_op_list)
#         rho_b = rho*n_b
#         rho_a = rho*sz
#         rho_c = rho*c
#         rho_d = rho*n_c
        
#         return rho_c.tr(),rho_a.tr(),rho_b.tr(),rho_d.tr()
    
    
#     if 'mapping' in kwargs:
               
#         H0 = H1 + H2w + H3 + H4 #+ H5
#         rho = steadystate(H0,c_op_list)
#         rho_c = rho*c
#         return rho_c.tr()
    
    if 'energies'in kwargs:
        H = H1 + H2 + H3
        
        return H.eigenenergies() #+ H4
    
    if 'energy'in kwargs:
        H = H1 + H2 + H3
        
        return H.eigenenergies(),s,alfa #+ H4
    
    if 'time'in kwargs:
        
        H0 = H1 + H2 + H3 #+ H4 
        H_args = {'H0': H0, 'c': c, 'cDag': c.dag() , 'A' : A , 'w': w, 's':s, 'sx':sx}

    #     rho = steadystate(H,c_op_list)

        T = 2 * pi / w

        U = propagator(Ht, T, c_op_list, H_args)

        rho_ss = propagator_steadystate(U)

        
        rho_b = rho_ss*n_b
        rho_a = rho_ss*sz
        rho_c = rho_ss*c + s*rho_ss*sx
        rho_d = rho_ss*n_c
        
        return rho_c.tr(),rho_a.tr(),rho_b.tr(),rho_d.tr()



In [149]:

    
N = 2
M = 2
P = 2
w_c = 5.0
w_nr = 3.5
g = 0.1
L = 0.001
Ej = 15
Ec = 0.223
w = 0
w_q_max = sqrt(8 * Ec * Ej) - Ec
w_q = 3.5
print(w_q_max)
d = 0.10
A = 0.00005# field aplitude
kwargs = {'num_cpus':26,'energy':1, 'Fred':1}

E,s,alfa= calc_spectrum_6(N,M,P, w_c,w_nr, w_q,L,g,A,w=0, **kwargs)









    



4.95000686255



In [150]:

    
E - E[0],s,alfa









    Out[150]:





(array([  0.        ,   3.49453058,   3.50000014,   5.00546914,
          6.99453044,   8.49806254,   8.50546928,  11.9980624 ]),
 -0.011764705882352941,
 -0.00014285714285714287)



In [172]:

    
N = 2
M = 2
P = 2
w_c = 5
w_nr = 3.5
g = 0.1
L = 0.001
Ej = 15
Ec = 0.223
w = 0
w_q_max = sqrt(8 * Ec * Ej) - Ec
print(w_q_max)
d = 0.10
A = 0.00005# field aplitude

# phi = linspace(0,pi/2,200)
x_i,x_f = 0.3,0.36
phi = pi*linspace(x_i,x_f,100)
x_vec =  sqrt( 8 * Ec * Ej* abs(cos(phi))*sqrt(1+(d*tan(phi))**2) )-Ec
# energies = array([calc_spectrum_6(N,M,P, w_c,w_nr, w_q/2,L,g,A,w,**kwargs)
#                   for w_q in x_vec])
kwargs = {'num_cpus':26,'energies':1, 'Fred':1}









    



4.95000686255



In [ ]:



In [173]:

    
# Energies


# variable to count the total number of tasks we need to do; used to create progress bar
task_count =len(x_vec)



# Check number of cpus to be used
if 'num_cpus' in kwargs:
    num_cpu = kwargs['num_cpus']
    if num_cpu == 1:
        print("1 CPU; Serial Simulation")
    else:
        print("Parallel Simulation with %d CPUs " % num_cpu)    
else:
    num_cpu = 1
    print("Serial Simulation")



## Program to run function in parallel: 

try:
    t_start = time.time() # start time simulation
    time_1 = []
    pool = mp.Pool(processes=num_cpu) #  create the initial pool to run the simulation   
#         manager = mp.Manager()
#         queue = manager.Queue()


#         _update_progress_bar(1)
#     task_args = a,z
    results = [pool.apply_async(calc_spectrum_6,(N,
                                                 M,
                                                 P,
                                                 w_c,
                                                 w_nr,
                                                 a1,
                                                 L,
                                                 g,
                                                 A,
                                                 w),kwargs
                                ,callback=None,error_callback=None) for a1 in x_vec]



    #####N,M,P, w_c,w_nr, w_q,L,g,A,w=0
    while True:
        incomplete_count = sum(1 for x in results if not x.ready())

        if incomplete_count == 0:
            print("[100.0%] of the simulations calculated, Estimated Remaining time: 0.0s", end="\r")
            print( "\nAll done! \nTotal time:%s"%datetime.timedelta(seconds=int(dif_time)))
            break

        else:

            p = float(task_count - incomplete_count) / task_count * 100 

            dif_time = (time.time() - t_start)    

#                 
            if p > 0:
                rem_time = (datetime.timedelta(seconds=int(dif_time*(100-p)/p)))

#                     rem_time_1 = (datetime.timedelta(seconds=int(dif_time/(task_count-incomplete_count))))
                time_1.append(float(dif_time/(task_count-  incomplete_count)))
#                     rem_time_1 = mean(time_1) *task_count
#                     rem_time_1 = (datetime.timedelta( seconds=int(mean(time_1) *task_count)))
                rem_time_1 = time.strftime("%Z - %Y/%m/%d, %H:%M:%S", time.localtime(t_start+mean(time_1) *task_count))
            else:
                rem_time = '?'
                rem_time_1 = 0


            print("[%4.1f%%] of the simulations calculated, Estimated Remaining time: %s, (%s)"
                  %(p,rem_time,rem_time_1) , end="\r")

            time.sleep(.25)


    while not all([ar.ready() for ar in results]):

        for ar in results:    
            ar.wait(timeout=0.1)

    pool.terminate()
    pool.join()

except KeyboardInterrupt as e:
    pool.terminate()
    pool.join()
    raise e




energies_temp = [ar.get() for ar in results]
energies = asarray(energies_temp)









    



Parallel Simulation with 26 CPUs 
[100.0%] of the simulations calculated, Estimated Remaining time: 0.0s
All done! 
Total time:0:00:00



In [174]:

    
# Plot 
fig, axes = subplots(1,1, figsize=(16,6))
x_inf = -1
x_sup = 10

for n in range(len(energies[0,:])):
    axes.plot(phi/pi, (energies[:,n]-energies[:,0]),'-',linewidth=2)
#     axes.plot(phi/pi, (energies[:,n]-energies[:,0])/2,'--')
    
    if n < 4:
        axes.text(x_i,energies[0,n]-energies[0,0],r'|%s>'%(n),fontsize=20)
    
axes.set_title('Full')
axes.set_ylim(x_inf, x_sup)
axes.set_xlabel(r'$\phi$', fontsize=18)
axes.set_ylabel(r'$E_n-E_0$', fontsize=18)
axes.hlines(w_nr,x_i,x_f,linestyles='dashed',linewidth=3,color='green')
axes.hlines(w_c,x_i,x_f,linestyles='dashed',linewidth=3)
axes.vlines(0.33,0,10,linestyles='dashed',linewidth=3)









    Out[174]:





<matplotlib.collections.LineCollection at 0x7fccca3a6ef0>



In [175]:

    
y_i,y_f = 4.9,5.1
y_vec = linspace(y_i,y_f,100) 
a , b = zip(*itertools.product(x_vec,y_vec))
kwargs = {'num_cpus':26,'time':1, 'Fred':1}



In [176]:

    
# Run Spectrum
# Create from the original vectors the new vector with the correct number copies
a , b = zip(*itertools.product(x_vec,y_vec))
# variable to count the total number of tasks we need to do; used to create progress bar
task_count =len(x_vec)*len(y_vec)

# Check number of cpus to be used
if 'num_cpus' in kwargs:
    num_cpu = kwargs['num_cpus']
    if num_cpu == 1:
        print("1 CPU; Serial Simulation")
    else:
        print("Parallel Simulation with %d CPUs " % num_cpu)    
else:
    num_cpu = 1
    print("Serial Simulation")



## Program to run function in parallel: 
t_start = time.time() # start time simulation
time_1 = []
try:
    pool = mp.Pool(processes=num_cpu) #  create the initial pool to run the simulation   
#         manager = mp.Manager()
#         queue = manager.Queue()


#         _update_progress_bar(1)
#     task_args = a,z
    results = [pool.apply_async(calc_spectrum_6,(N,
                                                 M,
                                                 P,
                                                 w_c,
                                                 w_nr,
                                                 a1,
                                                 L,
                                                 g,
                                                 A,
                                                 b1),kwargs
                                ,callback=None,error_callback=None) for a1,b1 in zip(a,b)]



        #####
    while True:
        incomplete_count = sum(1 for x in results if not x.ready())

        if incomplete_count == 0:
            print("[100.0%] of the simulations calculated, Estimated Remaining time: 0.0s", end="\r")
            print( "\nAll done! \nMean time:%f"%(dif_time/task_count))
            print( "\nTotal time:%s"%datetime.timedelta(seconds=int(dif_time)))
            break

        else:

            p = float(task_count - incomplete_count) / task_count * 100 

            dif_time = (time.time() - t_start)    

#                 
            if p > 0:
                rem_time = (datetime.timedelta(seconds=int(dif_time*(100-p)/p)))

#                     rem_time_1 = (datetime.timedelta(seconds=int(dif_time/(task_count-incomplete_count))))
                time_1.append(float(dif_time/(task_count -  incomplete_count)))
#                     rem_time_1 = mean(time_1) *task_count
#                     rem_time_1 = (datetime.timedelta( seconds=int(mean(time_1) *task_count)))
                rem_time_1 = time.strftime("%Z - %Y/%m/%d, %H:%M:%S", time.localtime(t_start+mean(time_1) *task_count))
            else:
                rem_time = '?'
                rem_time_1 = 0


            print("[%4.1f%%] of the simulations calculated, Estimated Remaining time: %s, (%s)"
                  %(p,rem_time,rem_time_1) , end="\r")

            time.sleep(.25)


    while not all([ar.ready() for ar in results]):

        for ar in results:    
            ar.wait(timeout=0.1)

    pool.terminate()
    pool.join()

except KeyboardInterrupt as e:
    pool.terminate()
    pool.join()
    raise e




results = [ar.get() for ar in results]









    



Parallel Simulation with 26 CPUs 
[100.0%] of the simulations calculated, Estimated Remaining time: 0.0s
All done! 
Mean time:0.294834

Total time:0:49:08



In [177]:

    
# Reshape Results
#results = qload('Two_Dispersive_Simulation')
results_1 = asarray(results)
# qsave(results,name='One_Dispersive_Simulation_200x300')
#qsave(results,name='Two_Dispersive_Simulation')
# qsave(results,name='ThirtytyVolts')

tr_c = reshape(results_1[:,0],(-1,len(y_vec+1)))
tr_a = reshape(results_1[:,1],(-1,len(y_vec+1)))
tr_b = reshape(results_1[:,2],(-1,len(y_vec+1)))
tr_d = reshape(results_1[:,3],(-1,len(y_vec+1)))



In [178]:

    
# Plot Graphics
fig, ax = subplots(4,1, figsize=(16,20))



im = ax[0].pcolor(phi/pi,y_vec,transpose(log10(abs(tr_c))))#,vmin=0, vmax=1)
fig.colorbar(im, ax=ax[0])
# ax[0,0].set_xlim(4.27,4.39)
# ax[0,0].set_ylim(P_i,P_f)
ax[0].set_ylabel(r'Probe Frequency (GHz)',fontsize=20)
# ax[0].set_xlabel(r'Qubit Tone Frequency (GHz)',fontsize=10)
ax[0].set_title(r'$Tr[\rho c]$',fontsize=20)


im = ax[1].pcolor(phi/pi,y_vec,transpose((abs(tr_a))))#,vmin=0, vmax=1)
fig.colorbar(im, ax=ax[1])
# ax[0,0].set_xlim(4.27,4.39)
# ax[0,0].set_ylim(P_i,P_f)
ax[1].set_ylabel(r'Probe Frequency (GHz)',fontsize=20)
# ax[1].set_xlabel(r'Qubit Tone Frequency (GHz)',fontsize=10)
ax[1].set_title(r'$Tr[\rho \sigma_z]$',fontsize=20)

im = ax[2].pcolor(phi/pi,y_vec,transpose(log10(abs(tr_b))))#,vmin=0, vmax=1)
fig.colorbar(im, ax=ax[2])
# ax[0,0].set_xlim(4.27,4.39)
# ax[0,0].set_ylim(P_i,P_f)
ax[2].set_ylabel(r'Probe Frequency (GHz)',fontsize=20)
# ax[2].set_xlabel(r'Qubit Tone Frequency (GHz)',fontsize=10)
ax[2].set_title(r'$Tr[\rho b^\dagger  b]$',fontsize=20)

im = ax[3].pcolor(phi/pi,y_vec,transpose((abs(tr_d))))#,vmin=0, vmax=1)
fig.colorbar(im, ax=ax[3])
# ax[0,0].set_xlim(4.27,4.39)
# ax[0,0].set_ylim(P_i,P_f)
ax[3].set_ylabel(r'Probe Frequency (GHz)',fontsize=20)
ax[3].set_xlabel(r'Flux ($\Phi_0$)',fontsize=20)
ax[3].set_title(r'$Tr[\rho c^\dagger c]$',fontsize=20)









    Out[178]:





<matplotlib.text.Text at 0x7fccca6cb4e0>



In [179]:

    
# Plot Graphic II
fig, axes = subplots(1,1, figsize=(16,10))
y_inf = y_i
y_sup = y_f
x_inf = x_i
x_sup = x_f

for n in range(len(energies[0,:])):
    axes.plot(phi/pi, (energies[:,n]-energies[:,0]),'-',linewidth=1)
    axes.plot(phi/pi, (energies[:,n]-energies[:,0])/2,'--')
    
#     if n < 4:
#         axes.text(.2,energies[0,n]-energies[0,0],r'|%s>'%(n),fontsize=20)
    
axes.set_title('Full')
axes.set_ylim(y_inf, y_sup)
axes.set_xlim(x_inf,x_sup)


axes.set_xlabel(r'$\phi$', fontsize=18)
axes.set_ylabel(r'Cavity Tone Frequency GHz', fontsize=18)
axes.hlines(w_nr,x_i,x_f,linestyles='dashed')
axes.hlines(w_c,x_i,x_f,linestyles='dashed')
# axes.vlines(0.245,0,10,linestyles='dashed',linewidth=3)
axes.vlines(0.33,0,10,linestyles='dashed')

im = axes.pcolor(phi/pi,y_vec,transpose(log10(abs(tr_c))))#axes.pcolor(phi/pi,y_vec,transpose((abs(tr))))#,vmin=0, vmax=1)
fig.colorbar(im, ax=axes)
# ax[0,0].set_xlim(4.27,4.39)
# ax[0,0].set_ylim(P_i,P_f)
# axes.set_ylabel(r'Qubit Tone Power(dBm)',fontsize=10)
# axes.set_xlabel(r'Qubit Tone Frequency (GHz)',fontsize=10)
# axes.set_title(r'$Tr[\rho\sigma_z]$',fontsize=20)









    Out[179]:





<matplotlib.colorbar.Colorbar at 0x7fccc9fc7278>

Scan Coupling



In [397]:

    
phi = 0.33 * pi
w_q = sqrt( 8 * Ec * Ej* abs(cos(phi))*sqrt(1+(d*tan(phi))**2) )-Ec


x_i,x_f = 0.0,0.005
x_vec= linspace(x_i,x_f,60)

y_i,y_f = 4.98,5.02
y_vec = linspace(y_i,y_f,30) 

a , b = zip(*itertools.product(x_vec,y_vec))
kwargs = {'num_cpus':26,'dispersive':1, 'Fred':1}



In [398]:

    
# Run Spectrum
# Create from the original vectors the new vector with the correct number copies
a , b = zip(*itertools.product(x_vec,y_vec))
# variable to count the total number of tasks we need to do; used to create progress bar
task_count =len(x_vec)*len(y_vec)

# Check number of cpus to be used
if 'num_cpus' in kwargs:
    num_cpu = kwargs['num_cpus']
    if num_cpu == 1:
        print("1 CPU; Serial Simulation")
    else:
        print("Parallel Simulation with %d CPUs " % num_cpu)    
else:
    num_cpu = 1
    print("Serial Simulation")



## Program to run function in parallel: 
t_start = time.time() # start time simulation
time_1 = []
try:
    pool = mp.Pool(processes=num_cpu) #  create the initial pool to run the simulation   
#         manager = mp.Manager()
#         queue = manager.Queue()


#         _update_progress_bar(1)
#     task_args = a,z
    results = [pool.apply_async(calc_spectrum_6,(N,
                                                 M,
                                                 P,
                                                 w_c,
                                                 w_nr,
                                                 w_q,
                                                 a1,
                                                 g,
                                                 A,
                                                 b1),kwargs
                                ,callback=None,error_callback=None) for a1,b1 in zip(a,b)]



        #####calc_spectrum_6(N,M,P, w_c,w_nr, w_q,L,g,A,w=0, **kwargs)
    while True:
        incomplete_count = sum(1 for x in results if not x.ready())

        if incomplete_count == 0:
            print("[100.0%] of the simulations calculated, Estimated Remaining time: 0.0s", end="\r")
            print( "\nAll done! \nMean time:%f"%(dif_time/task_count))
            print( "\nTotal time:%s"%datetime.timedelta(seconds=int(dif_time)))
            break

        else:

            p = float(task_count - incomplete_count) / task_count * 100 

            dif_time = (time.time() - t_start)    

#                 
            if p > 0:
                rem_time = (datetime.timedelta(seconds=int(dif_time*(100-p)/p)))

#                     rem_time_1 = (datetime.timedelta(seconds=int(dif_time/(task_count-incomplete_count))))
                time_1.append(float(dif_time/(task_count -  incomplete_count)))
#                     rem_time_1 = mean(time_1) *task_count
#                     rem_time_1 = (datetime.timedelta( seconds=int(mean(time_1) *task_count)))
                rem_time_1 = time.strftime("%Z - %Y/%m/%d, %H:%M:%S", time.localtime(t_start+mean(time_1) *task_count))
            else:
                rem_time = '?'
                rem_time_1 = 0


            print("[%4.1f%%] of the simulations calculated, Estimated Remaining time: %s, (%s)"
                  %(p,rem_time,rem_time_1) , end="\r")

            time.sleep(.25)


    while not all([ar.ready() for ar in results]):

        for ar in results:    
            ar.wait(timeout=0.1)

    pool.terminate()
    pool.join()

except KeyboardInterrupt as e:
    pool.terminate()
    pool.join()
    raise e




results = [ar.get() for ar in results]









    



Parallel Simulation with 26 CPUs 
[100.0%] of the simulations calculated, Estimated Remaining time: 0.0s
All done! 
Mean time:0.114948

Total time:0:03:26



In [399]:

    
#results = qload('Two_Dispersive_Simulation')
results_2 = asarray(results)
# qsave(results,name='One_Dispersive_Simulation_200x300')
#qsave(results,name='Two_Dispersive_Simulation')
# qsave(results,name='ThirtytyVolts')

tr_c = reshape(results_2[:,0],(-1,len(y_vec+1)))
tr_a = reshape(results_2[:,1],(-1,len(y_vec+1)))
tr_b = reshape(results_2[:,2],(-1,len(y_vec+1)))
tr_d = reshape(results_2[:,3],(-1,len(y_vec+1)))



In [400]:

    
# Plot Graphics
fig, ax = subplots(4,1, figsize=(16,20))



im = ax[0].pcolor(x_vec,y_vec,transpose(log10(abs(tr_c))))#,vmin=0, vmax=1)
fig.colorbar(im, ax=ax[0])
# ax[0,0].set_xlim(4.27,4.39)
# ax[0,0].set_ylim(P_i,P_f)
ax[0].set_ylabel(r'Probe Frequency (GHz)',fontsize=20)
# ax[0].set_xlabel(r'Qubit Tone Frequency (GHz)',fontsize=10)
ax[0].set_title(r'$Tr[\rho c]$',fontsize=20)


im = ax[1].pcolor(x_vec,y_vec,transpose((abs(tr_a))))#,vmin=0, vmax=1)
fig.colorbar(im, ax=ax[1])
# ax[0,0].set_xlim(4.27,4.39)
# ax[0,0].set_ylim(P_i,P_f)
ax[1].set_ylabel(r'Probe Frequency (GHz)',fontsize=20)
# ax[1].set_xlabel(r'Qubit Tone Frequency (GHz)',fontsize=10)
ax[1].set_title(r'$Tr[\rho \sigma_z]$',fontsize=20)

im = ax[2].pcolor(x_vec,y_vec,transpose(log10(abs(tr_b))))#,vmin=0, vmax=1)
fig.colorbar(im, ax=ax[2])
# ax[0,0].set_xlim(4.27,4.39)
# ax[0,0].set_ylim(P_i,P_f)
ax[2].set_ylabel(r'Probe Frequency (GHz)',fontsize=20)
# ax[2].set_xlabel(r'Qubit Tone Frequency (GHz)',fontsize=10)
ax[2].set_title(r'$Tr[\rho b^\dagger  b]$',fontsize=20)

im = ax[3].pcolor(x_vec,y_vec,transpose((abs(tr_d))))#,vmin=0, vmax=1)
fig.colorbar(im, ax=ax[3])
# ax[0,0].set_xlim(4.27,4.39)
# ax[0,0].set_ylim(P_i,P_f)
ax[3].set_ylabel(r'Probe Frequency (GHz)',fontsize=20)
ax[3].set_xlabel(r'$\lambda$ GHz',fontsize=20)
ax[3].set_title(r'$Tr[\rho c^\dagger c]$',fontsize=20)









    Out[400]:





<matplotlib.text.Text at 0x7f10d711ea20>